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Abstract 

Our  previous  PTTI  paper  demonstrated  that  a  two-state  Kalman  filter  involving  the 
parameters  frequency  and  frequency  drift,  when  properly  implemented,  produces  frequencies 
for  hydrogen  masers  that  generate  a  stable  mean  timescale,  at  least  for  cases  involving 
simulated  and  low-noise  real  data.  The  current  paper  extends  the  investigation  to  the  noisier 
data  actually  employed  in  the  USNO  operational  timescale.  Results  for  about  500  days  of 
postprocessed  data  for  12  masers  and  51  cesium  clocks  indicate  comparable  frequency 
stabilities  to  those  obtained  with  the  currently  operational  timescale  algorithm. 


1  INTRODUCTION 

Mean  timescales  based  on  Kalman  filter  modelling  of  individual  clocks  have  been  investigated  by  others 
in  order  to  provide  an  ensemble  time  that  is  more  uniform  than  the  time  kept  by  any  one  of  the  constituent 
clocks  and  whose  estimated  states  are  optimal  in  the  minimum  squared  error  sense.  Tryon  and  Jones  [1], 
Barnes  and  Allan  [2],  Stein  [3],  and  Greenhall  [4]  employed  as  state  parameters  the  phase,  frequency,  and 
frequency  drift  of  an  ensemble  of  cesium  clocks  relative  to  a  master  clock.  A  problem  with  this  type  of 
filter  is  that,  since  time  is  unobservable,  elements  of  the  covariance  matrix  grow  without  bound  [5]. 
Weiss  and  Weissert  [5]  utilized  a  one-state  filter  to  determine  the  frequency  of  cesium  clocks  used  in  their 
AT2  timescale.  Only  Tryon  and  Jones  and  Weiss  and  Weissert  tested  their  algorithms  on  real  data. 

Barnes  and  Allan  [2]  proposed  a  Kalman  filter  that  utilized  only  frequency  and  drift,  pointing  out  that 
only  these  two  parameters  are  physically  meaningful  for  a  frequency  standard,  and  indeed  showed  by 
simulations  that  such  a  two-parameter  filter  performed  better  than  the  three-parameter  one.  Since  the 
phase  (time)  of  a  clock  is  integrated  from  the  frequency,  it  adds  no  information,  and  solution  for  a  phase 
parameter  is  unnecessary. 

Consequently,  a  reasonable  model  for  both  our  cesium-beam  frequency  standards  and  our  hydrogen 
masers  (which  generally  have  inherent  frequency  drift)  is: 

f  (t  +  At )  =  fit')  +  At  •  r/(()  +  £(t  +  At) 
c/((  +  A()=  d(t)  +  ij{t  +  A() 


511 


Report  Documentation  Page 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 


1.  REPORT  DATE 

DEC  2002 


2.  REPORT  TYPE 


3.  DATES  COVERED 

00-00-2002  to  00-00-2002 


5a.  CONTRACT  NUMBER 


5b.  GRANT  NUMBER 


5c.  PROGRAM  ELEMENT  NUMBER 


5d.  PROJECT  NUMBER 


5e.  TASK  NUMBER 


5f.  WORK  UNIT  NUMBER 


4.  TITLE  AND  SUBTITLE 

Kalman  Filter  Characterization  of  Cesium  Clocks  and  Hydrogen  Masers 


6.  AUTHOR(S) 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES)  8.  PERFORMING  ORGANIZATION 

U.S.  Naval  Observatory ,3450  Massachusetts  Avenue  report  number 

NW, Washington, DC, 20392-5420 

9.  SPONSORING/MONITORING  AGENCY  NAME(S )  AND  ADDRESS(ES )  10.  SPONSOR/MONITOR' S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

34th  Annual  Precise  Time  and  Time  Interval  (PTTI)  Planning  Meeting,  3-5  December  2002,  Reston,  VA 

14.  ABSTRACT 

Our  previous  PTTI  paper  demonstrated  that  a  two-state  Kalman  filter  involving  the  parameters  frequency 
and  frequency  drift,  when  properly  implemented,  produces  frequencies  for  hydrogen  masers  that  generate 
a  stable  mean  timescale,  at  least  for  cases  involving  simulated  and  low-noise  real  data.  The  current  paper 
extends  the  investigation  to  the  noisier  data  actually  employed  in  the  USNO  operational  timescale.  Results 
for  about  500  days  of  postprocessed  data  for  12  masers  and  51  cesium  clocks  indicate  comparable 
frequency  stabilities  to  those  obtained  with  the  currently  operational  timescale  algorithm. 

15.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

18.  NUMBER 

19a.  NAME  OF 

ABSTRACT 

OF  PAGES 

RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Same  as 
Report  (SAR) 

16 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


34th  Annual  Precise  Time  and  Time  Interval  (PTTI)  Meeting 


where  /  (t)  is  a  clock's  frequency  at  time  t  =  1,2,  d  (t)  is  its  frequency  drift;  s  (t)  and  ?/  (t)  are 
independent  random  variables  with  zero  mean  and  normal  distributions  uncorrelated  in  time  (zero 
autocorrelation),  i.e.  white  noise  processes;  and  At  is  the  time  step.  The  latter  is  1  hour  in  our  case,  which 
is  in  the  white  FM  noise  regime  of  cesiums  and  in  the  flicker  FM  regime  of  masers.  For  the  masers, 
measurement  system  noise  dominates  clock  noise  for  sampling  times  up  to  a  few  hours  and  is 
characterized  by  white  phase  noise.  A  Kalman  filter  that  assumes  the  presence  of  only  white  FM  and 
random-walk  FM  noises  cannot  exactly  model  1  -hour  data  from  a  maser.  Flowever,  in  our  previous  paper 
[61,  it  was  shown  that,  for  an  ensemble  of  hydrogen  masers  connected  to  a  low-noise  (Timing  Solutions 
Corporation  or  TSC)  measurement  system,  such  a  filter  with  the  parameters  of  frequency  and  frequency 
drift  is  quite  capable  of  producing  a  mean  timescale  with  a  frequency  stability  that  is  not  only  viable,  but 
at  sampling  times  up  to  several  days,  actually  superior  to  that  produced  by  a  conventional  least-squares 
algorithm  such  as  the  currently  operational  at  USNO  [7].  Among  other  issues,  this  paper  addresses 
whether  a  viable  timescale  can  be  similarly  generated  for  the  noisier  data  from  our  currently  operational 
Data  Acquisition  System. 

In  order  to  minimize  process  noise,  one  must  choose  a  reference  that  is  stable  as  possible.  The  Mean 
timescale,  at  least  as  derived  from  the  clocks  under  consideration,  cannot  so  serve  because  it  is  not  yet 
available  at  the  first  step  of  a  recursion.  One  could  refer  all  the  frequencies  and  drifts  to  one  of  the  clocks 
or  to  a  Mean  predicted  from  previous  data.  Flowever,  it  would  be  easier  and  more  error-  and  correlation- 
free  to  refer  one  type  of  clocks,  say  the  masers,  to  the  Mean  of  the  other  type  of  clocks  (cesiums);  this  is 
our  objective.  In  USNO’s  operational  timescale,  the  maser  rates  and  drifts  are  calibrated  against  the 
cesium  Mean  anyway,  which  is  always  available. 


2  THE  KALMAN  FILTER 

The  state  transition  equations  are: 
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X(t  +  At)  =  0(0  X(t)  +  W  (t  +  At) 

where  the  subscripts  denote  clocks  and  O  (t)  is  the  state  transition  matrix. 

Let  Z;  (t)  be  the  phase  of  clock  i  at  time  t  relative  to  the  same  reference.  The  observation  equations  are: 
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or: 

Z(t,M)  =  H(t)  X(t)  +  V(t) 

where  Z  ( t,At )  is  the  vector  of  measurements,  H  (t)  is  the  observation  matrix  (time-dependent  in  that  the 
number  of  weighted  clocks  with  available  data  may  vary),  and  V  (t)  is  the  vector  of  measurement  errors  V; 
(0- 
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where  S  is  the  noise  spectral  density  (sometimes  called  variance  rate)  of  random-walk  FM  and  s  is 

the  noise  spectral  density  of  random-run  FM  (random  walk  of  frequency  drift).  The  effect  of  these  errors 
on  the  system  for  any  interval  At  is: 

Q(At)  =  f  <D(f)  ^  0  T(t)  dt 
J  dt 

0 
0 
0 
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where  O7  is  the  transpose  of  matrix  O. 
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For  uncorrelated  parameters,  the  covark 
small  At,  is  such  that: 
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If  the  measurement  errors  are  uncorrelated,  the  measurement  error  vector  V  (t  +  At)  has  the  covariance 
matrix: 


R(t) 


0  Jv  2  0 

v2 

o  o  '• 


where  sv2  is  the  noise  spectral  density  of  random  walk  of  phase,  or  white  FM. 

The  assumption  of  zero  autocorrelation  does  require  prefdtering  for  such  outliers  as  erroneous  phase 
measurements  shared  by  consecutive  first  differences.  The  real  data  in  this  and  the  previous  paper  were 
subjected  to  robust  prefiltering. 

Rather  than  noise  spectral  densities,  let  us  instead  solve  for  the  variances  oy  =  sv  At,  ae  =  se  At,  and  o,f  = 
s,,2At  as  our  noise  parameters.  The  Hadamard  variance,  Oh,  when  the  noise  sources  consist  only  of  white 
noise  FM,  random-walk  FM,  and  random-run  FM,  can  be  expressed  as  [8]: 

(T/=(Tvy  +  (1/6)  cr/  r  +  (11/120)  o2  r3  (i) 

where  x  is  the  sampling  time.  Our  noise  parameters  may  be  solved  by  least  squares  from  this  equation  of 
condition.  Hadamard  variance  has  the  beneficial  property  of  being  insensitive  to  frequency  drift. 

If  the  reference  were  a  mean  timescale,  improved  values  could  be  obtained  from  fits  of  the  individual 
clocks  and  recomputed  occasionally  at  times  when  the  performance  of  a  clock  appears  to  have  changed 
significantly. 

The  first  step  in  the  version  of  the  Kalman  recursion  in  which  no  knowledge  of  the  process  noise  is 
required  is  to  calculate  the  inverse  of  the  transition  covariance: 

P~\t)  =  [p-{t)]A+HT{t)RA{t)H{t) 

where  the  transition  covariance  matrix  predicted  for  the  next  step  P~(t)  (whose  superscript  denotes 
“predicted”)  can  be  assumed  initially  to  contain  huge  covariances. 

Second,  we  compute  the  Kalman  gain  K: 

K  (t)  =  P  (t)  HT (t)  R~l(t) 

where  P  ( t )  is  the  transition  covariance  matrix  of  the  current  step.  Next,  we  update  the  parameter  estimate 
thusly  (where  A  denotes  “estimated”): 

x(t)  =  r  (/)+£(/)  [Z(t)-H(t)r(t)} 

Then  we  make  the  following  predictions  for  the  next  step 
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r(t+ ao  =  0(0  x(t) 

p-{t+M )  =  0(0  ?(0  or(0  +  Qit) 

and  so  on.  Each  application  of  the  recursion  yields  an  estimate  of  the  system  state  that  is  a  function  of  the 
elapsed  time  since  the  last  fdter  update,  which  can  occur  any  time,  i.e.  At  is  not  necessarily  constant,  and 
the  data  need  not  be  equally  spaced.  Note  that  there  are  two  matrix  inversions  instead  of  one  as  in  the 
more  common  formulation  of  the  Kalman  filter.  We  saw  in  our  previous  paper  that  this  is  of  no 
consequence  [6], 

If  prefilter  analysis  of  the  data  indicates  the  presence  of  a  significant  frequency  step,  the  covariances  in 
P  (t)  could  be  reset  to  huge  values  in  order  to  reinitialize  the  filter.  As  formulated  here,  the  Kalman  filter 
regains  the  state  information  within  a  few  steps,  and  the  data  are  appropriately  downweighted  until  that 
has  occurred.  This  resetting  of  the  covariances  affects  only  the  clock  involved,  not  the  entire  timescale, 
which  is  generated  from  averaging  the  Kalman  frequencies  across  the  entire  ensemble,  as  explained  in  the 
next  section. 


3  THE  KALMAN  MEAN 

The  frequencies  and  drifts  produced  by  this  two-state  Kalman  filter  could  be  averaged  and  integrated  to 
generate  a  Kalman-filter-based  mean  timescale.  The  relative  sizes  of  the  clock  covariances  P~(t)  would 
be  inversely  related  to  the  individual  clock  weights.  Robustness  would  require  the  imposition  of  an  upper 
limit  on  any  one  clock's  relative  weight.  If  the  filter  and  model  operate  correctly,  the  clock  variances 
should  slowly  converge  to  steady-state  values. 

While  there  is  some  degradation  of  the  frequency  stability  of  the  timescale  as  a  result  of  the  capping  of 
each  clock’s  relative  weight,  this  degradation  is  rarely  significant  in  practice  unless  the  ensemble  is  small 
or  very  inhomogeneous,  neither  of  which  is  anticipated  to  be  the  case  upon  implementation.  In  any  case, 
robustness  must  be  purchased  at  some  price,  and  this  price  will  have  to  be  decided  upon  during 
implementation. 

If,  in  forming  the  mean  timescale,  each  clock’s  frequency  is  weighted  according  to  the  inverse  of  its 
variance,  the  timescale  would  be  biased  by  the  “clock-ensemble  effect.”  The  prediction  error  is  always 
too  small  because  the  frequency  of  a  clock  is  correlated  with  the  frequency  of  the  ensemble,  since  the 
ensemble  includes  a  contribution  from  each  clock.  The  weights,  which  are  proportional  to  the  inverse  of 
the  variances,  are  therefore  systematically  too  large,  causing  a  positive  feedback  that  increasingly  biases 
the  timescale  toward  the  ensemble’s  best  clock  [9,10].  This  can  be  avoided  if  the  weight  of  each  clock  is 
based  on  its  stability  relative  to  a  Mean  of  which  it  is  not  a  part,  either  by:  (1)  using  another  Mean 
entirely,  say  one  based  on  cesium  clocks  rather  than  masers;  or  (2)  computing  a  Mean  consisting  of  the 
rest  of  the  clocks  in  the  ensemble  and  referring  the  clock  in  question  to  that  Mean.  (2)  can  be 
accomplished  by  using  the  relation: 

au2  =  a2  /  (1  -  w) 

where  o2  is  the  uncorrected  variance,  au2  is  the  corrected  (unbiased)  variance,  and  weight  is  the  relative 
weight,  assumed  to  be  constant  with  time  [7].  In  practice,  an  upper  limit  would  also  have  to  be  placed  on 
any  one  clock’s  weight. 
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Parameter  and  error  estimates  would  be  available  from  the  filter  in  near  real  time,  but  practical 
implementation  would  require  robust  detection  and  rejection  of  outliers  (measurements  whose  errors  are 
unlikely  to  have  originated  in  the  clock),  as  well  as  prompt  recognition  of  time  and  frequency  steps. 
Kalman  filters  provide  both  a  forecast  of  the  next  data  point  and  an  estimate  of  its  uncertainty,  so  it  is 
possible  to  implement  robust  outlier  detection,  though  one  must  allow  for  the  possibility  that  discrepant 
points  are  indications  of  a  step  in  time  or  frequency.  Stein  proposed  an  altered  gain  function  to  smoothly 
deweight  outliers  and  an  adaptive  filter  to  recognize  time  and  frequency  steps  [3].  Data  step  detection 
will  be  addressed  in  Section  6. 

Initially  the  filter's  knowledge  of  the  parameters  will  be  crude,  but  in  postprocessing  one  can  make  use  of 
“future  data”  by  running  forward-  and  backward-moving  filters  through  the  data  and  averaging  the 
corresponding  parameters  weighted  by  the  number  of  time  steps  contributing  to  each  using  weights  based 
on  filter  variances  (“smoothing”)  [11].  In  order  to  avoid  using  the  same  datum  twice,  the  state  at  any  one 
time  from  one  of  the  filters  has  to  be  averaged  with  the  predicted  state  at  that  time  from  the  other  filter. 
This  kind  of  processing  of  time  data  has  only  been  done,  to  our  knowledge,  by  NIST  a  posteriori  for  their 
AT2  timescale  [12].  This  could  be  done  in  near  real  time  every  time  step  for  some  interval  of  the  latest 
data,  though  retroactive  recomputation  would  be  involved. 


4  THE  MASERS 

The  two-state  Kalman  filter  described  above,  including  the  forward-and-backward  averaging  and  a 
correction  for  the  clock-ensemble  effect,  was  applied  to  539  days  of  hourly  data  selected  from  those  of  12 
Sigma  Tau/Datum/Symmetricom  masers  at  USNO.  The  data  were  produced  by  our  recently  renovated 
Data  Acquisition  System  (DAS),  a  measurement  system  involving  universal  time-interval  counters  and 
multiplexed  switches.  The  dominant  noise  of  this  system  is  a  phase  noise  of  ~35  ps. 

The  data  were  referenced  to  the  Mean  of  the  same  12  masers  as  computed  by  the  current  USNO  timescale 
algorithm  [7],  which  removed  constant  frequencies  and  drifts  from  each  steady-state  segment  of  every 
maser  by  modelling  it  against  the  cesium  Mean  computed  by  the  same  algorithm.  In  practice,  the 
incestuous  nature  of  the  ensemble  and  its  reference  can  be  avoided  by  referencing  the  masers  to  the  Mean 
predicted  by  the  filter  and  modelling  the  masers  against  a  Kalman  cesium  Mean.  But  in  this  paper,  an 
independent  reference  is  not  necessary,  since  our  objectives  are  to:  (1)  test  the  stability  and  lack  of 
systematic  error  of  the  Kalman  Mean  relative  to  the  operational  Mean  for  DAS  data  of  both  masers  and 
cesiums,  rather  than  the  lower-noise  maser  data  investigated  in  our  previous  paper  [6];  (2)  determine  the 
time  constant  of  the  filter  in  modelling  simulated  shocks  to  the  clocks,  which  will  involve  only 
comparisons  between  shocked  and  unshocked  data;  and  (3)  search  for  an  optimal  clock  weighting 
scheme,  whose  relative  values  should  not  change  significantly  between  two  strongly  correlated  references 
(i.e.  an  independent  cesium  Mean  and  a  cesium-calibrated  Mean). 

The  Kalman  filter  yielded  very  satisfactory  results  for  the  DAS  data  of  each  maser.  For  example,  Figure 
1  compares  the  Kalman  frequencies  with  the  hourly  frequencies  of  phase  for  maser  NAVI 2.  The  537 
days  of  data  shown  would  have  to  have  been  least-squares  fit  by  more  than  five  linear  segments  in  order 
to  be  detrended  for  frequency  and  drift  with  our  current  USNO  algorithm  [7].  Deciding  on  the  placement 
of  such  segments  requires  labor  and  some  subjective  judgment  that  the  Kalman  filter  permits  us  to  avoid. 

The  Kalman  frequencies  were  not  very  sensitive  to  choice  of  process  noise  parameters  (tuning  of  the 
filter),  i.e.  the  fit  of  each  clock’s  entire  Hadamard  variance  curve  by  Equation  (1)  gave  sufficiently 
accurate  noise  parameters  in  spite  of  the  occasional  small  changes  in  the  rate  and  drift  parameters 
encountered  here.  The  noise  parameters  for  the  masers  had  the  following  ranges: 
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crv2  =2.0'10“30  f>2.3'10"29 
c/  =  5.8'10"33  o5.5'10“30 
=  1.2  10“37  <-»2.3'l(T34 

The  phase  error  of  the  DAS  measurement  system  is  ~35  ps,  so  for  an  hourly  first  difference,  the 
observational  noise  would  be  (V2  •  35  ps/hr)2  =  1.9  •  10‘2S. 

5  THE  KALMAN  MASER  MEAN 

Generation  of  a  mean  timescale  by  the  simple  integration  of  the  average  instantaneous  Kalman  frequency 
estimates  for  the  clock  ensemble  yields  a  pretematurally  stable  result  reflecting  only  the  evolution  of 
one’s  clock  models.  Connection  with  reality  must  be  maintained  through  use  of  the  basic  timescale 
equation  [7,13]: 


*?  =  zt~A,  +  w o  wo  -  wo  +  ft( o  a*  -  } 

1  2 

where  z  is  the  phase  of  the  reference  minus  the  Mean  of  the  ensemble;  wt  (?)  is  the  weight  of  clock  i; 
X.  (?)  is  the  phase  of  the  reference  minus  clock  i;  f  (?)  is  the  frequency  and  dt  (?)  is  the  frequency  drift 

of  clock  i,  respectively,  relative  to  the  Mean,  all  at  time  t;  and  Aris  the  time  step.  In  our  case,  A.I  =  1 
hour.  Both  the  Kalman  Mean  and  its  reference,  the  current  USNO  Mean,  were  calculated  with  this 
equation;  the  only  difference  is  that  f  (?)  and  d  t  (?)  are  the  state  parameters  of  the  filter  and  the  products 

of  least-squares  fits  to  first  differences,  respectively.  Unless  noted  otherwise,  the  individual  clock 
weights  are  unity  unless  a  data  point  is  rejected  in  the  prefiltering. 

One  of  the  computational  advantages  to  averaging  clock  frequencies  (as  above)  instead  of  clock  phases 
(as  in  timescale  algorithms  that  use  another  form  of  the  basic  equation)  is  the  obviation  of  the  time  steps 
in  the  mean  timescale  when  weighted  clocks  are  added  or  removed  from  the  ensemble.  In  the  case  of  the 
current  USNO  timescale  algorithm,  the  clock’s  rate  relative  to  the  Mean  is  determined  prior  to  its  addition 
to  the  ensemble  in  order  that  this  rate  can  be  removed  by  the  basic  equation,  preserving  the  (assumed) 
constant  frequency  of  the  ensemble.  In  the  case  of  the  Kalman  filter  under  test,  however,  the  clock 
frequency  is  determined  relative  to  the  current  USNO  maser  Mean,  rather  than  the  Kalman  Mean  being 
computed,  so  the  rate  to  be  removed  is  that  relative  to  the  USNO  maser  Mean,  to  which  the  Kalman  Mean 
will  then  be  referenced.  This  will  not  be  the  case  in  an  actual  implementation,  however.  In  that  case,  the 
cesiums  will  be  referenced  to  the  Kalman  cesium  Mean  being  computed,  and  the  masers  will  be  calibrated 
by  and  referenced  to  the  Kalman  cesium  Mean. 

Averaging  the  12  masers  together  over  the  time  span  MJD  52048-52587  (though  some  of  the  masers 
came  in  as  late  as  MJD  52347),  we  obtained  the  resulting  Kalman  Mean,  relative  to  the  USNO  maser 
Mean,  is  shown  in  Figure  2.  The  two  Means  do  not  differ  more  than  1.4  ns  from  one  another.  Since  both 
Means  have  been  generated  from  the  same  raw  data  and  the  clocks  in  each  Mean  were  ultimately 
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calibrated  against  the  USNO  cesium  Mean,  they  are  of  very  similar  long-term  stability,  and  their  phases 
should  random-walk  away  from  one  another,  which  indeed  appears  to  be  the  case. 

In  order  to  determine  the  absolute  frequency  stabilities  of  the  two  sets  of  maser  Means  (current  and 
Kalman),  we  divided  the  masers  into  three  sub-ensembles  of  four  masers  each  and  performed  a  three- 
cornered-hat  analysis  on  each  set  using  the  method  of  Torcaso  et  al.  [14],  which  corrects  for  sub-ensemble 
intercorrelations.  The  results  are  given  in  Figure  3.  The  Kalman  timescale’s  stabilities  are  at  least  as 
great  as  those  of  the  current  timescale’s.  Apparently  the  significant  measurement  (phase)  noise  at  the 
shortest  sampling  times  and  the  flicker  FM  noise  of  the  masers  are  not  impediments  to  the  use  of  this 
filter  on  masers,  though  the  separation  between  the  two  sets  of  stability  curves  at  the  mid-term  sampling 
times  in  Figure  3  indicates  a  different  apportioning  of  the  process  and  measuring  errors  by  the  two 
timescale  algorithms. 

In  order  to  use  data  efficiently  and  combine  clocks  properly  in  a  Kalman  mean,  we  should  test  possible 
clock  weighting  schemes.  The  criterion  for  such  a  scheme  would  be  maximum  stability  in  the  sampling 
time  region  of  most  interest  for  the  steering  of  our  Master  Clock,  which  ranges  from  1  day  (the  time 
constant  of  our  frequency  steering)  to  25  days  (the  time  constant  of  our  phase  steering)  [15].  It  might 
seem  logical  to  use  as  clock  weights  the  inverse  frequency  variances  estimated  by  the  Kalman  filter. 
However,  it  is  not  possible  to  investigate  weighting  schemes  using  these  data,  because  they  were 
calibrated  with  a  Mean  based  on  the  same  (equally  weighted)  clocks,  so  any  stability  difference  between 
the  Kalman  and  current  Means  due  to  the  weighting  would  only  indicate  that  their  respective  weighting 
schemes  were  different.  As  we  will  see  in  Section  7  for  the  cesiums,  weights  based  on  frequency 
stabilities  for  sampling  times  shorter  than  the  time  constant  of  the  filter  (found  to  be  1 1  days  in  Section  6) 
yield  better  results  than  those  based  on  longer  sampling  times.  But  the  case  for  the  masers  is  complicated 
by  the  fact  that  masers  are  already  beyond  the  white  FM  noise  regime  at  our  shortest  sampling  time  (1 
hour).  Indeed,  whiteness  tests  with  the  MATLAB  software  package  indicate  that:  (1)  the  frequency 
residuals  of  the  individual  masers  are  significantly  autocorrelated  at  sampling  times  shorter  than  1 6  hours 
due  to  each  clock’s  prevailing  flicker  noise;  and  (2)  the  frequency  residuals  from  the  Kalman  maser  Mean 
are  significantly  autocorrelated  at  sampling  times  shorter  than  16  days,  probably  due  to  the  fact  that  the 
masers  are  being  calibrated  against  an  incestuous  Mean  (one  computed  from  the  same  data  processed  by 
another  timescale  algorithm)  and  to  the  1 1  -day  time  constant  of  the  filter. 


6  KALMAN  FILTER  TIME  CONSTANT 

Before  moving  on  to  the  noisier  data  of  the  cesiums,  let  us  use  the  maser  data  to  determine  the  time 
constant  of  this  Kalman  filter  when  operating  on  DAS  data.  This  time  constant  is  important  because  of 
the  effects  of  process  noise,  such  as  that  due  to  environmental  shocks  to  the  clocks.  Significant  shocks  on 
longer  timescales  than  this  time  constant  can  be  properly  handled  by  the  filter  in  its  normal  operation,  but 
such  shocks  that  take  place  on  shorter  timescales  would  require  special  handling,  e.g.  data  rejection  if  the 
effect  is  transient,  or  clock  downweighting  and  covariance  reinitialization  if  the  effect  appears  to  be 
permanent,  at  least  until  the  filter  can  again  satisfactorily  model  the  clock’s  new  parameters.  Failure  to 
take  appropriate  action  could  result  in  the  destabilization  of  the  Master  Clock  if  it  were  steered  to  the 
Kalman-filter-based  Mean. 

In  view  of  this  eventuality,  before  this  filter  can  be  implemented  in  real  time,  a  test  will  need  to  be 
designed  that  will  recognize  frequency  excursions  in  the  data  that  are  too  rapid  to  be  followed  closely  by 
the  filter.  Perhaps  the  frequency  innovations  produced  by  a  Kalman  filter  tuned  for  white  phase  and  white 
FM  noises  (and,  hence,  having  a  much  shorter  time  constant  than  that  of  our  proposed  filter),  with  any 
significant  frequency  divergences  taken  as  indicating  a  recent  frequency  step  by  the  clock.  Another 
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possibility,  at  least  for  the  masers,  is  to  monitor  their  frequencies  on  a  system  with  a  lower  measuring 
noise  than  the  DAS,  e.g.  the  TSC  system  used  in  our  previous  study  [6];  cesiums  are  not  profitably 
studied  by  such  a  system  because  of  their  greater  process  noise  exceeds  the  measuring  noise  of  even  the 
DAS. 

Prior  to  developing  such  a  test,  we  must  determine  the  time  constant  of  the  current  filter.  This  can  be 
done  by  inserting  an  artificial  step  in  the  first- difference  data  and  observing  how  long  the  filter 
frequencies  require  to  return  to  steady  state.  Ramp  shocks  of  various  sizes  and  over  intervals  from  1  hour 
to  1  day  were  inserted  in  the  data  for  each  of  the  12  masers.  A  sample  of  the  results  is  shown  in  Figure  4 
for  four  of  the  masers.  The  average  recovery  time  was  found  to  be  11  days,  which  we  will  take  as  the 
time  constant  of  this  filter. 

Consequently,  the  filter  requires  a  rate-change  detector  with  a  time  constant  much  shorter  than  1 1  days. 


7  THE  CESIUMS  AND  THE  KALMAN  CESIUM  MEAN 

Having  tested  our  filter  on  the  stabler  data  of  the  masers,  let  us  now  proceed  to  the  noisier  data  of  the 
cesium  clocks,  which  must  ultimately  provide  the  systematic  calibration  of  all  the  clocks.  As  stated 
earlier,  however,  in  this  test  we  will  calibrate  the  frequencies  of  the  cesiums  against  those  of  the  masers, 
whose  frequencies  were  in  turn  calibrated  with  those  of  these  very  cesiums,  so  there  will  be  incest  not 
present  in  the  envisioned  implementation.  Doing  so  for  498  days  of  data  for  51  cesium  clocks,  we 
obtained  the  following  ranges  for  the  process  noise  parameters: 

<7V2  =  3.7  •  10“28  <-a  1 .0  •  10“27 
a/  =  2.7'10"32  <-»1.6'10"29 
^  =  4.M0'36  o  9.4-10'33 

The  noise  model  assumed  by  our  filter  is  valid  for  the  cesiums  because  the  DAS  measurement  phase 
noise,  even  had  it  been  a  problem  for  the  masers,  is  swamped  by  the  white  FM  process  noise  of  the 
cesiums. 

Assuming  equal  clock  weights  for  the  moment,  we  computed  a  Kalman  cesium  Mean  whose  differences 
from  the  current  USNO  cesium  mean  are  given  in  Figure  5  and  whose  frequency  stability  relative  to  the 
USNO  maser  Mean  is  depicted  in  Figure  6.  Figure  5  shows  no  departure  of  the  two  Means  greater  than  8 
ns,  about  what  the  means  would  be  expected  to  depart  from  one  another  as  the  result  of  random  walk  FM 
[16].  In  Figure  6,  we  see  that  the  two  Means  are  of  comparable  stability  on  the  long  term.  The  small 
decrease  in  stability  on  the  short  term  is  likely  due  to  the  more  realistic  (as  opposed  to  optimistic) 
modelling  of  the  white  noise  of  the  individual  cesiums  compared  to  the  least-squares  fit  of  the  current 
USNO  algorithm.  In  Section  5  we  also  noted  effects  from  the  different  noise  modelling  for  the  masers. 

Unlike  the  masers,  the  cesiums  still  display  white  FM  noise  for  sampling  times  as  long  as  our  1-hour 
measurement  rate,  so  it  might  be  possible  to  increase  the  stability  of  the  Kalman  cesium  Mean  by 
weighting  the  clocks  unequally.  We  have  not  been  able  to  accomplish  this  with  the  current  USNO 
algorithm,  partly  because  the  clocks  are  so  similar  in  frequency  stability  and  partly  because  any  clock 
displaying  a  change  in  its  parameters  is  immediately  deweighted  until  it  can  be  remodelled,  a  procedure 
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that  tends  to  make  the  performance  of  all  weighted  clocks  appear  similar.  Still,  perhaps  the  more  realistic 
treatment  of  the  process  noises  by  the  Kalman  fdter  may  make  differentiation  between  the  clocks 
possible. 

The  Kalman  frequencies  possess  the  dominant  clock  noise,  so  it  is  logical  to  try  to  weight  the  clocks 
according  to  the  inverse  frequency  variances  estimated  every  time  step  by  the  filter.  Weighting  the  clocks 
thusly,  however,  only  worsened  the  stability  (see  Figure  7),  probably  because  the  frequency  variances  are 
mostly  a  measure  of  the  noise  prevailing  around  sampling  times  similar  to  the  time  constant  of  the  filter 
(11  days),  while  the  relative  stabilities  of  the  individual  clocks  evidently  change  more  rapidly. 

Perhaps,  then,  it  might  be  possible  to  increase  the  stability  of  the  Kalman  cesium  Mean  by  weighting  the 
clocks  according  to  the  inverse  Allan  variances  for  sampling  times  much  shorter  than  1 1  days. 
Consequently,  the  cesium  data  were  weighted  using  Allan  variances  for  sampling  times  of  1  hour,  6 
hours,  12  hours,  and  1  day.  Indeed,  the  shorter  the  sampling  time,  the  better  the  stability  became  (see 
Figure  7),  but  the  stability  did  not  quite  attain  that  of  our  equally  weighted  results.  Consequently,  it  is  not 
likely  that  one  can  improve  the  stability  of  this  Kalman-filter-based  Mean  with  Allan-variance  weights. 


CONCLUSION 

We  have  shown  that: 

•  the  DAS  measurement  data  for  the  masers,  the  same  used  in  the  operational  USNO  maser  mean 
timescale,  can  also  be  successfully  utilized  by  our  Kalman  filter  to  generate  a  Mean  that  is  both  as 
stable  as  the  operational  maser  Mean  and  as  that  produced  from  the  less  noisy  TSC  data  studied 
previously  [6],  at  least  over  sampling  times  including  a  few  weeks,  which  is  necessary  if  such  a 
timescale  were  to  replace  the  current  one  for  purposes  of  steering  a  Master  Clock. 

•  DAS  measurement  data  for  the  cesium  clocks,  the  same  used  in  the  operational  USNO  cesium  mean 
timescale,  can  also  be  successfully  utilized  by  our  Kalman  filter  to  generate  a  cesium  Mean  about  as 
stable  as  the  operational  cesium  Mean,  within  the  limitations  of  the  incestuous  frequency  calibration. 
This  capability  is  necessary  if  such  a  timescale  were  used  to  replace  the  current  one  for  purposes  of 
calibrating  the  frequencies  and  drifts  of  the  masers  in  our  ensemble. 

•  The  time  constant  of  our  filter/data  combination  is  1 1  days,  meaning  that  significant  shocks  to  the 
clocks  occurring  over  intervals  shorter  than  1 1  days  must  be  isolated  and  their  effects  eliminated  from 
the  Mean  using  a  method  with  a  more  alacritous  step  response  than  that  of  our  filter. 

No  better  weighting  scheme  was  found  for  the  cesiums  than  that  of  equality  among  each. 

Before  our  final  objective,  real-time  implementation  of  this  mean  timescale  algorithm,  can  be  attempted, 

we  must: 

•  test  and  implement  a  method  capable  of  identifying  and  allowing  for  state  changes  too  rapid  for  this 
proposed  filter  to  handle; 

•  test  whether  a  free-running  Mean  predicted  by  this  algorithm  can  serve  as  a  sufficiently  stable 
reference,  rather  than  the  operational  USNO  maser  Mean  used  in  this  paper; 
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•  decide  whether  averaging  the  state  parameters  from  the  forward-  and  backward-moving  filters  is 
preferable  to  using  those  of  the  forward-moving  filter,  given  that  the  smoothing  by  the  former  is 
purchased  at  the  cost  of  phase  steps  in  the  real-time  Mean. 

The  advantages  of  replacing  our  current  timescale  algorithm  with  a  Kalman-filter-based  one  would  be  a 
reduction  in  labor  and  an  increase  in  the  mathematical  rigor  of  the  Mean  computation. 
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Figure  1.  A  comparison  between  hourly  frequencies  and  Kalman  frequencies  for  maser  NAV12  relative 
to  the  current  USNO  maser  Mean. 


Figure  2.  The  phase  difference  between  the  Kalman  maser  Mean  and  the  USNO  maser  Mean.  The  clocks 
are  equally  weighted  in  both. 


523 


34,h  Annual  Precise  Time  and  Time  Interval  (PTTI)  Meeting 


LOG  TAU  (sec) 


CURRENT 

ALGORITHM 

ENSEMBLE  #1 

CURRENT 

ALGORITHM 

ENSEMBLE  #2 

CURRENT 

ALGORITHM 

ENSEMBLE  #3 

KALMAN 

ALGORITHM 

ENSEMBLE  #1 

KALMAN 

ALGORITHM 

ENSEMBLE  #2 

KALMAN 

ALGORITHM 

ENSEMBLE  #3 


Figure  3.  Frequency  stability  of  maser  Means  for  three  sub-ensembles  as  generated  by  the  Kalman  filter 
and  the  current  USNO  timescale  algorithm.  The  error  bars  were  computed  using  the  method  of  Ekstrom 
and  Koppang  [17]. 
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Figure  4.  Kalman  filter  frequencies  before  and  after  insertion  of  a  simulated  1  ns/day  step  at  MJD 
523 12.0  for  masers  NAV2  (red),  NAV4  (green),  NAV8  (blue),  and  NAVI  1  (orange). 
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Figure  5.  The  phase  difference  between  the  Kalman  cesium  Mean  and  the  USNO  maser  Mean.  The 
clocks  are  equally  weighted  in  both. 
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Figure  6.  The  frequency  stability  vs.  sampling  time  x  for  the  Kalman  and  USNO  cesium  Means  relative 
to  the  USNO  maser  Mean. 
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Figure  7.  The  frequency  stability  vs.  sampling  time  x  for  Kalman  cesium  Means,  relative  to  the  USNO 
maser  Mean,  computed  for  equal  clock  weights  and  clock  weights  based  on  Allan  variances  for  sampling 
times  of  1  hour,  6  hours,  12  hours,  and  1  day. 
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